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We present a new determination of the Bs leptonic decay constant from lattice QCD simulations 
that use gluon configurations from MILC and a highly improved discretization of the relativistic 
quark action for both valence quarks. Our result, fss =0.225(4) GeV, is almost three times more 
accurate than previous determinations. We analyze the dependence of the decay constant on the 
heavy quark's mass and obtain the first empirical evidence for the leading l/y/ruh dependence 
predicted by Heavy Quark Effective Theory (HQET). As a check, we use our analysis technique 
to calculate the ttiBs — mass difference. Our result agrees with experiment to within errors 

of llMeV (better than 2%). We discuss how to extend our analysis to other quantities in Bs and 
B physics, making 2%-precision possible for the first time. 



PACS numbers: 11.15.Ha,12.38.Aw,12.38.Gc 

Lattice simulations of QCD have become essential 
for high-precision experimental studies of 5-meson de- 
cays — studies that test our understanding and the lim- 
itations of the standard model of weak, electromagnetic 
and strong interactions, and also determine fundamental 
parameters, like the CKM matrix, in that model. Ac- 
curate theoretical calculations of QCD contributions to 
meson masses, decay constants, mixing amplitudes, and 
semileptonic form factors are critical for this program, 
and lattice simulation is the main tool for providing these 
calculations. A major complication for the lattice sim- 
ulations has been the large mass of the h quark, which 
has necessitated the use of non-relativistic effective field 
theories like NRQCD to describe h dynamics in the sim- 
ulations. The need for effective field theories has made it 
difficult to achieve better than 5-10% precision for many 
important quantities. 

Recently we overcame the analogous problem for 
c quarks by introducing a highly improved discretization 
of the relativistic quark action that gives accurate re- 
sults even on quite coarse lattices: the Highly Improved 
Staggered-Quark (HISQ) discretization [l]. With this 
formalism, c quarks are analyzed in the same way as li, 
(i, and s quarks, which greatly reduces the uncertainties 
in QCD simulations of D physics [l-Q- More recently 
we showed that the HISQ action can be pushed to much 
higher masses — indeed, very close to the h mass — us- 
ing new lattices, from the MILC collaboration, with the 
smallest lattice spacing available today (a = 0.044 fm). 
This allowed us to extract a value for the 6's MS mass 
that was accurate to better than 1%. Here we extend 
that work in a new analysis of the Bs meson's leptonic 
decay constant /^s, which produces the most accurate 
theoretical value to date. 

We also compute the mass difference ttib^ — ^77^/2, as 



an additional test of our analysis method. This differ- 
ence is particularly sensitive to QCD dynamics because 
the leading (and uninteresting) dependence on the heavy 
quark's mass mostly cancels in the difference. 

It would be quite expensive to extend our new analy- 
sis directly to 5-meson quantities, because of the added 
costs associated with very light valence and sea quarks. 
This is unnecessary, however, because heavy-quark ef- 
fective field theories like NRQCD already give accurate 
results for ratios of Bs to B quantities, like fss/fB- This 
is because the largest systematic errors from these effec- 
tive theories, due to operator matching, cancel in such 
ratios. The ratio of the decay constants, for example, is 
known to ±2% from effective field theories [8]. So the 
combination of accurate Bs quantities from HISQ simu- 
lations, as discussed in this paper, with Bs/ B ratios from 
NRQCD or other effective field theories provides a potent 
new approach to high-precision h physics generally. Note 
that no operator matching is required in our relativistic 
analysis of fs^ because of the exact chiral symmetry of 
the HISQ formalism in the massless limit. 

In our simulations for this paper, we computed decay 
constants and masses for non-physical Hs mesons com- 
posed of an s quark, and heavy quarks h with various 
masses rrih ranging from below the c mass to just be- 
low the h mass. This data allows us to extrapolate to 
the h mass, where ttiHs =^Bs- We repeated our analy- 
sis for five different lattice spacings, allowing us also to 
extrapolate our results to zero lattice spacing. 

The gluon configuration sets we used are from the 
MILC collaboration [9] and are described in Table HI 
Our simulation results for the decay constants and meson 
masses, for various values of the h mass, are presented 
in Table |TT1 We also give results for the mass of the 
pseudoscalar hh meson, m^^^ , and for the mass of the rjs 
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FIG. 1. The leptonic decay constant for pseudoscalar Ks mesons ifs, plotted on the left versus the Hs mass as the h- 
quark's mass is varied. The solid line and gray band show our best-fit estimates for the decay constants extrapolated to zero 
lattice spacing. Best-fit results (dashed lines) and simulation data are also shown for five different lattice spacings, with results 
for smaller lattice spacings extending to higher masses (since we restrict amu < !)• The simulation data points have been 
corrected for small mistunings of the s quark's mass. On the right the same simulation data and fits are plotted for ^^ttiHs Jhs 
versus I/ttihs • 



TABLE I. Parameter sets used to generate the 3-flavor gluon 
configurations analyzed in this paper. The lattice spacing 
is specified in terms of the static-quark potential parameter 
ri =0.3133(23) fm [10|]; values for n/a are from [9|]. The bare 
quark masses are for the ASQTAD formalism and uq is the 
fourth root of the plaquette. The spatial (L) and temporal 
(T) lengths of the lattices are also listed, as are the number 
of gluon configurations (A^cf) and the number of time sources 
(Ats) per configuration used in each case. 



Set 


n/a 


auomou/d 


auomos 


uo 


L/a 


T/a 


A^ef X Nts 


1 


2.152(5) 


0.0097 


0.0484 


0.860 


16 


48 


631 X 2 


2 


2.618(3) 


0.01 


0.05 


0.868 


20 


64 


595 X 2 


3 


3.699(3) 


0.0062 


0.031 


0.878 


28 


96 


566 X 4 


4 


5.296(7) 


0.0036 


0.018 


0.888 


48 


144 


201 X 2 


5 


7.115(20) 


0.0028 


0.014 


0.895 


64 


192 


208 X 2 



meson. The rjs is an unphysical pseudoscalar ss meson 
whose valence quarks are not allowed to annihilate; we 
use its mass to tune the bare mass of the s quark: simu- 
lations show that its mass is rrin^^phys = 0.6858(40) GeV 
when the s mass is correctly tuned |lOj . 



TABLE IL Simulation results for each of the five configura- 
tion sets (Table [J) and several values of the heavy-quark's 
mass TTih- The s-quark's mass nis is tuned to be close to its 
physical value. Results are given for: the leptonic decay con- 
stant and mass tuhs of the pseudoscalar hs meson, and 
masses of the pseudoscalar hh and ss mesons, rrir^^ and rrir^^ 
respectively. 



aMr, 



aruh aMHs 



1 0.061 
0.061 



0.5049(4 
0.5049(4 



2 0.0492 
0.0492 
0.0492 



0.4144(2 
0.4144(2 
0.4144(2 



0.0337 
0.0337 
0.0337 
0.0337 



0.2941(1 
0.2941(1 
0.2941(1 
0.2941(1 



0.0228 
0.0228 
0.0228 
0.0228 



0.2062(2 
0.2062(2^ 
0.2062(2^ 
0.2062(2 



0.0165 
0.0165 
0.0165 
0.0165 



0.1548(1 
0.1548(1 
0.1548(1 
0.1548(1 



0.66 
0.81 



1.3108(6 
1.4665(8 



0.44 
0.63 
0.85 



0.9850(4 
1.2007(5 
1.4289(8 



0.3 
0.413 
0.7 
0.85 



0.7085(2 
0.8472(2 
1.1660(4 
1.3190(5 



0.273 
0.564 
0.705 
0.85 



0.5935(2 
0.9313(5 
1.0811(8 
1.228(1) 



0.195 
0.5 
0.7 
0.85 



0.4427(3 
0.8038(8 
1.017(1) 
1.168(2) 



0.1913(7) 
0.197(1) 



1.9202(2) 
2.1938(2) 



0.1500(5) 
0.1559(7) 
0.161(1) 



1.4240(1) 
1.8085(1) 
2.2193(1) 



0.1054(2) 
0.1084(2) 
0.1112(5) 
0.1123(6) 



1.03141(8 
1.28057(7" 
1.86536(5 
2.14981(5 



0.0750(3) 
0.0754(6) 
0.0747(8) 
0.074(1) 



0.8994(1) 
1.52542(6 
1.80845(6 
2.08753(6 



0.0555(3) 
0.055(1) 
0.053(2) 
0.052(2) 



0.67113(6 
1.34477(8 
1.75189(7^ 
2.04296(7 



We expect some statistical correlation between results 
from the same configuration set but with different h- 
quark masses. We have not measured these, but we have 
verified that our results are insensitive (at the level of 
=bcr/4) to such correlations. We introduce a 50% correla- 
tion for our fits, which increases our final error estimates 
slightly. 

Our strategy for extracting fs^ is first to fit our simu- 



lation results for fn^ to the HQET-inspired formula [11 1 

fH,{a,mH,,mrjJ = 



/ / X \ -2//3o Nm-1 / 



+ C,(m^^ -^,,phys). (1) 

where /3o = ll — 2n//3 = 9 in our simulations [l2[, ay is 
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the QCD couphng in the V scheme [H, [l3|, and constant 
6 = — 0.5 from HQET. Here we use rrtH^ as a proxy for 
the /i-quark mass since its value for b quarks is known 
from experiment. The last term in Eq. ([T]) corrects for 
tuning errors in the 5-quark mass; we determined that 
Cs = 0.06(1) by repeating our calculations with slightly 
different s masses [l^. We parameterize 




(2) 

with Nm = Na = 4:. This expansion is in powers of quark 
masses and the QCD scale parameter Aqcd ^ 0.5 GeV 
divided by the ultraviolet cutoff for the lattice theory: 
Auv ^Tv/a. The fit parameters are the coefficients Cijki 
for each of which we use a prior of ± 1.5, which is 
conservative [15]. 

Our data for five different lattice spacings and a wide 
range of masses are presented with our fit results 
in Fig [H The reach in grows as the lattice spac- 
ing decreases (since we restrict arrih < 1), and deviations 
from the continuum curve get smaller. The fit is excel- 
lent, with a per degree of freedom of 0.36 while fitting 
all 17 measurements. The small results from our con- 
servative priors (we get excellent fits and smaller errors 
with priors that are half the width). 

Having determined the parameters in Eq. ([T]), the sec- 
ond step in our analysis is to set Mh^ =Mb^^ a = 0, and 
rrirj^ =^r7s,phys 1^ that formula to obtain our final value 
for /b,, 

/b^ = 0.225(4) GeV, (3) 

which agrees well with the previous best NRQCD re- 
sult of 0.231(15) GeV [16] but is almost four times more 
accurate. Our result also agrees with the recent result 
of 0.232(10) GeV from the ETM cohaboration, although 
that analysis includes only two of the three light quarks 
in the quark sea [17]. 

The error budget for our result is given in Table [ml and 
shows that the dominant errors come from statistical un- 
certainties in the simulations, the -^'^Bs extrapo- 
lation, the a^^O extrapolation, and uncertainties in the 
scale-setting parameter ri. Our analysis of f^^ i^i [6] in- 
dicates that finite volume errors, errors due to mistuned 
sea-quark masses, errors from the lack of electromagnetic 
corrections, and errors due to lack of c quarks in the sea 
are all significantly less than 1%, and so negligible com- 
pared with our other uncertainties. Our final result is 
also insensitive to the detailed form of the fit function; 
for example, doubling the number of terms has negligible 
effect (0.03cr) on the errors and value. 

We have also included in Fig. [T] (right) a plot of 
y/rriHs for different values of mH^' This shows that 
there are large non- leading terms in fn^i beyond the 
leading \j ^mn^ behavior predicted by HQET. Our sim- 
ulation nevertheless provides evidence for the leading 



TABLE ni. Dominant sources of uncertainty in our determi- 
nations of the Bs decay constant and the Bs — r\h mass dif- 
ference. Contributions are shown from the extrapolations in 
vfiHs 5 and vfis , as well as statistical errors in the simulation 
data and errors associated with the scale-setting parameter 
r\. Other errors are negligible. 







rriBs - ^r7b/2 


Monte Carlo statistics 


1.30% 


1.49% 


VTiHs ruBs extrapolation 


0.81 


0.05 


ri uncertainty 


0.74 


0.33 


^ extrapolation 


0.63 


0.76 


mr^g ^r7s,phys extrapolation 


0.13 


0.18 


n/a uncertainties 


0.12 


0.17 


Total 


1.82% 


1.73% 



term. Treating exponent b in Eq. ([T]) as a fit parameter, 
rather than setting it equal to —0.5, we find a best-fit 
value of b= —0.51(13), in excellent agreement with the 
HQET prediction. This is the first empirical evidence for 
this behavior. 

Our analysis also yields a value for f^^ , which agrees 
with [6]. It is also clear from Fig. [T] (left) that fn^ peaks 
between fjj^ and /b^, and that /b^ is smaller — we find: 

=0.906(14). (4) 

HQET suggests a ratio less than one, but previous lattice 
QCD results have been ambiguous about this point. 

To check our /bs analysis technique, we adapted the 
same technique to compute the mass difference [18] 

A = mH,-mr^j2, (5) 

using, as inputs, the masses m^^^ computed in our sim- 
ulations for pseudoscalar hh mesons made of our heavy 
quark. Our values for A come from the results in Ta- 
ble im We fit them to 

A{a,mH,,mr^J 

(6) 

where our simulations indicate that ds = 0.18(1) [l9| . 
and Di{a) has an expansion similar to that for Ci{a) 
(Eq. ([2])), with the same priors. 

We show the best fit to our simulation data in Fig. [2l 
Again, the fit is excellent, which a P^r degree of free- 
dom of 0.13 while fitting all 17 measurements. Extrapo- 
lating to the b mass, we obtain our best fit value for the 
mass splitting, 

rriB^ - mr^j2 = 0.658(11) GeV, (7) 

which agrees well with experiment: experiment gives 
0.671(2) GeV which becomes 0.666(4) after removing 



4 



CD 



0.70 


.1 


1 1 


1 1 


0.65 








0.60 








0.55 
0.50 
0.45 
0.40 


/ mo, 
1 


1 1 


i 

1 




2 


3 4 
m^, (GeV) 


5 



FIG. 2. The Hs — rju/'^ mass difference plotted versus mi/^ 
as the h quark's mass is varied. The sohd hne and gray band 
show our best-fit estimates for the mass differences extrapo- 
lated to zero lattice spacing. Best-fit results (dashed lines) 
and simulation data are also shown for five different lattice 
spacings, with results for smaller lattice spacings extending 
to larger masses (since we require arrih < !)• The simulation 
data points have been corrected for small mistunings in the 
s quark's mass. Data points (in black) at rriDs and ttiBs are 
the experimental values after correcting for small effects from 
electromagnet ism, rfh annihilation, and c quarks in the sea, 
none of which are included in the simulation. 



corrections from electromagnetism, 775 annihilation, and 
c quarks in the sea (not included in our simulations) 0, 
[20]. Our fit also gives a value for — it also 

agrees well with experiment [6]. 

In this paper, we have shown how to use a highly 
improved discretization of the relativistic quark action 
to make accurate calculations for mesons containing 
h quarks. Our result for the Bs decay constant, /b^ = 
0.225(4) GeV, agrees with other determinations from lat- 
tice QCD but is almost three times more accurate than 
the most precise previous result. The reliability of our 
extrapolations is underscored both by our previous de- 
termination of the 6-quark's MS mass, which agrees with 
other determinations to within errors of less than ±1%, 



and by our calculation here of the ttib^ 



difference, which agrees with experiment to within errors 
of ±llMeV or less than 2%. Our analysis of the decay 
constant gives the most extensive information to date on 
the heavy-quark mass dependence of the decay constant, 
and provides the first empirical evidence for the leading 
1/ -s/rrih dependence predicted by HQET. Further results 
on the Be mesons, and the decay constant fr^^ will be 
presented elsewhere [2l| . 

Our analysis has important implications for future lat- 
tice simulations of B physics. Other Bs quantities, like 
semileptonic form factors, can be analyzed in the same 
way, bringing few-percent precision within reach. Simi- 
lar precision for B quantities is possible by combining Bs 



calculations like these with precise calculations of Bs/B 
ratios using (very efficient) non-relativistic effective field 
theories for 6-quark dynamics. 
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